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Reconstruction of Graph Signals through 
Percolation from Seeding Nodes 

Santiago Segarra, Antonio G. Marques, Geert Leus, and Alejandro Ribeiro 


Abstract —New schemes to recover signals defined in the 
nodes of a graph are proposed. Our focus is on reconstructing 
bandlimited graph signals, which are signals that admit a sparse 
representation in a frequency domain related to the structure 
of the graph. Most existing formulations focus on estimating 
an unknown graph signal by observing its value on a subset 
of nodes. By contrast, in this paper, we study the problem of 
reconstructing a known graph signal using as input a graph 
signal that is non-zero only for a small subset of nodes (seeding 
nodes). The sparse signal is then percolated (interpolated) across 
the graph using a graph filter. Graph filters are a general¬ 
ization of classical time-invariant systems and represent linear 
transformations that can be implemented distributedly across 
the nodes of the graph. Three setups are investigated. In the 
first one, a single simultaneous injection takes place on several 
nodes in the graph. In the second one, successive value injections 
take place on a single node. The third one is a generalization 
where multiple nodes Inject multiple signal values. For noise¬ 
less settings, conditions under which perfect reconstruction is 
feasible are given, and the corresponding schemes to recover 
the desired signal are specified. Scenarios leading to imperfect 
reconstruction, either due to insufficient or noisy signal value 
injections, are also analyzed. Moreover, connections with classical 
interpolation in the time domain are discussed. The last part 
of the paper presents numerical experiments that illustrate the 
results developed through synthetic graph signals and two real- 
world signal reconstruction problems: influencing opinions in a 
social network and inducing a desired brain state in humans. 

Index Terms —Graph signal processing. Signal reconstruction. 
Interpolation, Percolation, Graph-shift operator, Bandlimited 
graph signals 


I. Introduction 

Sampling and reconstruction of bandlimited signals are 
cornerstone problems in classical signal processing. The emer¬ 
gence of new fields of knowledge such as network science 
and big data is generating a pressing need to extend the 
results existing for classical time-varying signals to signals 
defined on graphs |^-Q. This not only entails modifying the 
existing algorithms, but also gaining intuition on the concepts 
that are preserved and lost when a signal is defined not 
in the time grid, but in a more general graph domain. In 
the context of reconstruction of graph signals, two different 
approaches have been developed. A first approach is related 
to the interpolation of bandlimited signals, which consists in 
inferring unobserved values by leveraging the fact that the 
signal lives in a low-dimensional space |6)-|8J. Although 
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most interpolation approaches are centralized, iterative Q 
and distributed o) interpolation schemes have also 

been developed. A different approach towards graph signal 
reconstruction is graph regularization GD.IID where a notion 
of smoothness is assumed on the signal and the unobserved 
values are estimated based on this notion. Both approaches co¬ 
incide in that they estimate a graph signal from the observation 
of a subset of the signal values. By contrast, the approach in 
this paper is to preserve the two-step methodology used when 
recovering bandlimited time-varying signals, which consists in 
the generation of a sparse signal followed by the application of 
a low-pass filter to reconstruct the missing entries, and extend 
it to the more general graph domain. 

To be more specific, we study the reconstruction of ban¬ 
dlimited graph signals through the application of low-pass 
graph filters to sparse seeding signals. Graph filters are the 
generalization of classical time-invariant systems when the 
signals are defined on a general graph as opposed to the 
classical time domain Q. Seeding signals are graph signals 
attaining nonzero values on a subset of the nodes in the graph, 
called seeding nodes. To describe our approach more precisely, 
let y stand for the target graph signal we want to recover. Our 
goal is to design a graph filter H and a sparse signal x such 
that y can be obtained upon applying H to x. The design 
is accomplished in two steps. In the first step, we design the 
filter H leveraging the bandlimitedness of y, to eliminate the 
frequencies not present in y. Then, we use the H designed 
in the first step and the specific value of y to design the 
signal X. The challenge is that x cannot be chosen freely 
but is rather the output of a seeding phase where only a few 
seeding nodes inject values. The seeding phase requires more 
elaboration than its counterpart for time-varying signals, not 
only because graph signals are less regular, but also because 
it will be shown that in general the seeding values cannot 
coincide with those of the signal to recover. For a rigorous 
problem definition see Section [lT^ and Figure[T]in Section [HI] 
Since graph filters act on graph signals through the successive 
application of local operators, the output of a graph filter 
can be viewed as the outcome of a diffusion or percolation 
process. Applications include the generation of an opinion 
profile in a social network m by influencing a few agents 
(Section VII-B i and the synthesis of brain signals Hg by 
exciting a few neural regions (Section VII-C| l. Other potential 
applications for signal reconstruction via local interactions 
include molecular communications in nanonetworks eg, (13 
and wireless sensor networks (H. 

The paper investigates three different reconstruction 
schemes, each of them associated with a different seeding 
phase. In Section III the seeding phase consists of a unique 
seeding signal with several nonzero values, which coincides 
with the intermediate signal x. By contrast, in Section IV the 
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seeding phase consists of several seeding signals injected by 
a single seeding node. At each instant, the signal is perco¬ 
lated (diffused) within one-hop neighborhoods. The support 
of X depends on the duration of the seeding phase and the 
connectivity of the seeding node. Finally, in Section |V] we 
consider a more general scheme which merges the two earlier 
approaches. In this scheme, the seeding phase consists of 
several time instants and, in each of them, multiple nodes are 
allowed to inject a signal. The schemes will be referred to as 
multiple node-single time (MN-ST), single node-multiple time 
(SN-MT) and multiple node-multiple time (MN-MT) seeding, 
respectively. For the three of them, we state conditions on the 
underlying graph and the seeding nodes to guarantee perfect 
reconstruction of any bandlimited signal. We also show that, 
in general, if the interpolator takes the form of a graph filter, 
the seeding values cannot coincide with those of the signal to 
interpolate. Furthermore, we discuss how additional seeding 
values can be used to reduce the complexity of the graph 
filter needed for perfect recovery and draw connections with 


classical interpolation of time-varying signals. In Section VI 


we study the reconstruction performance in imperfect settings, 
either because the seeding values are insufficient in number 
or corrupted by noise. In Section VII we run numerical 


experiments to illustrate signal reconstruction in noiseless and 
noisy scenarios using both synthetic and real-world graphsj^ 


II. Bandlimited graph signals and graph filters 

Let Q denote a directed graph with a set of N nodes or 
vertices N and a set of links £, such that if node i is connected 
to j, then {i,j) G E. The (incoming) neighborhood of i is 
defined as the set of nodes Mi = {j \ {j,i) G £} connected 
to i. For any given graph we define the adjacency matrix A 
as a sparse N x N matrix with nonzero elements Aji if and 
only if {i,j) G £. The value of Aji captures the strength of 
the connection from i to j. The focus of this paper is not on 
analyzing Q, but a graph signal defined on the set of nodes 
M. Formally, such a signal can be represented as a vector 
X = [xi, G where the i-th component represents 

the value of the signal at node i or, alternatively, as a function 

f : M U, defined on the vertices of the graph. 

The graph ^ is endowed with a graph-shift operator S Q, 
m- The shift S is a N x N matrix whose entry Sji can be 

nonzero only if i = j or if {i, j) G £ . The sparsity pattern of 

the matrix S captures the local structure of Q, but we make no 
specific assumptions on the values of the nonzero entries of S. 
Choices for S are the adjacency matrix of the graph Q, P). 
its Laplacian 0, and their respective generalizations p^ . 
The intuition behind S is to represent a linear transformation 
that can be computed locally at the nodes of the graph. More 
rigorously, if y is defined as y = Sx, then node i can compute 
Pi provided that it has access to the value of Xj at j G Mi- We 
assume henceforth that S is diagonalizable, so that there exists 
a N X N matrix V and a N x N diagonal matrix A that can 
be used to decompose S as S = VAV~^. In particular, S is 
diagonalizable when it is normal, i.e., it satisfies SS^ = S^S 

* Notation: is the ith iV x 1 canonical basis vector (all entries of e; are 

zero except the ith one, which is one); Ejc '■= [ei,^..,e^] is a tall matrix 
collecting the K first canonical basis vectors while Ek '■= [oic-i-li ■••iGiv] 
collects the last N — K canonical basis vectors; 0 and 1 are, respectively, 
the all-zero and all-one matrices (when not clear from the context, a subscript 
indicating the dimensions will be used). 


where denotes the conjugate transpose of S. In that case, 
we have that V is unitary, which implies = V^, and 

leads to the decomposition S = VAV^. 

We are interested in cases where the graph-shift operator S 
plays a role in explaining the graph signal x. More specifically, 
cases where x can be expressed as a linear combination of a 
subset of the columns of V = [vi,..., Vjy], or, equivalently, 
where the vector x = V“^x is sparse |21|. In this context, 
vectors are interpreted as the graph frequency basis, Xi as 
the corresponding signal frequency coefficients, and x as a 
bandlimited graph signal. We assume that the set of active 
frequencies are known and, without loss of generality, that 
those are the first K ones associated with the eigenvalues of 
largest magnitude J2T), 1^. Under this assumption, if we 
denote by := [xi,..., xk]"'" a AT x 1 vector collecting the 
coefficients corresponding to those frequencies, it holds that 
X is a AT-bandlimited signal if 


X = [x^,0,... ,0]^, X = Vx := VifXif, (1) 

where we have defined the tall matrix := [vi,...,vj^] 
containing the first K eigenvectors of the shift operator S. 


A. Graph filters 


Graph filters H : -G are linear graph-signal 

operators of the form H := hiS^; i.e., polynomials (of 

degree A — 1) of the graph-shift operator Q. A particularity of 
graph filters is that they can be implemented locally, e.g., with 
A — 1 exchanges of information among neighbors. This is true 
because the application of S on a signal x can be computed 
through local interactions. 

The graph filter H can also be written as H = 
V(V“^. The diagonal matrix H := hi-S} 

can then be viewed as the frecjuency resgonse of H and it 
can be alternatively written as H = diag(h), where vector h 
is a vector that contains the frequency response of the filter. 
Let \i denote the Ath eigenvalue of S and define the N x L 
Vandermonde matrix 


/ 1 Ai ... Af-i \ 
\ 1 \n ■ • ■ A ^“^ / 


( 2 ) 


Upon defining the vector containing the coefficients of the 
filter as h := [hg,... it holds that h = and 

therefore 


H =J2fSg hiS‘ =Vdiag(^h)V-i=Vdiag(h)V-A (3) 


This implies that if y is defined as y = Hx, its frequency 
representation y satisfies 


y = diag(’®'h)x. 


(4) 


Within this context, a low-pass graph filter of bandwidth K is 
one where the frequency response h := fFh is given by 

h= [h^,0,...,0]^, (5) 

where h.K contains the frequency response for the first K 
frequencies. Notice that when the low-pass filter in 0 is 
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applied to an arbitrary signal x, the output signal is K- 
bandlimited as described in Q- An alternative expression to 
define a graph hlter is | [23j 

L-l 

H = ao n(S-ail), (6) 

1=1 


in y. This phase has duration L—1, which is the order of the 
hlter H. 

The goal of this paper is to design and H such 

that z = y. In Sections |I^|IV] and|V]we present this design for 
three different seeding schemes, where we impose additional 
restrictions on the structure and the number of seeding signals. 


which also gives rise to a polynomial on S of degree L — 1. A 
specihc advantage of the representation in ^ is that it provides 
a straightforward way to design low-pass hlters via successive 
annihilation of graph frequencies. In particular, if we hx o/ = 
Afe for some eigenvalue of S then the hlter H will eliminate 
the frequency basis v^, i.e., the eigenvector associated with A^. 
For future reference, we denote by D the number of distinct 
eigenvalues in 

Remark 1 (Discrete-time signals) To establish connections 
with classical time-varying signals, we dehne the directed 
cycle graph Qdc, with node set Af = { 1 , 2,. .., iV} and edge 
set £dc — {(*,modAr(f) + l)}i^i, where modAr(f) denotes 
the remainder obtained after dividing i by N. Its adjacency 
and Laplacian matrices are denoted, respectively, as Adc and 
Ijdc'=i^A.dc- Discrete-time periodic signals can be thought 
as graph signals on the directed cycle Qdc- Setting the shift 
operator either to S = Adc or S = hdc gives rise to the 
Fourier basis F. More formally, the right eigenvectors of S 
satisfy V = F, with Fij := exp(-|-j27r(f — l)(j — l)/N)/\/N 
where j := y/—l. Selecting S = Adc has the additional 
advantage of satisfying = exp(—j27r(f — l)/iV), i.e., the 
eigenvalues of the shift operator correspond to the classical 
discrete frequencies. Interpretations for the eigenvalues of hdc 
also exist Q. The frequency representation x of a graph signal 
X is given by x = V~^x whereas the frequency response of 
a filter with coefficients h is given by h = ’®'h. For general 
graphs, matrices and need not be related. However, 
for the case of Qdc, if S^c = Adc, then 'F = y/NF^ and 
. Thus, the Fourier transforms for signals and hlter 
coefficients are equivalent up to a constant for time-varying 
signals but this is not true for general graph signals. 


B. Signal reconstruction using graph filters 

Our objective is to reconstruct a specihc AT-bandlimited 
signal y by applying a graph hlter H to a signal x, where 
X is the result of a seeding procedure. More specihcally, the 
reconstruction scheme proceeds in two phases; 

• Seeding phase. The input to this phase is a set of t 

graph signals denominated seeding signals. 

These signals percolate through the graph following the 
dynamics 

x(*) = -b x(-i) = 0. (7) 

The output of this phase is set as x := 

• Filtering phase. The graph signal x is used as input to a 
low-pass graph hlter H, generating the output z := Hx. 

The purpose of the seeding phase, which has duration r, is to 
inject into the graph the information needed to interpolate the 
signal y. The hltering phase further propagates the informa¬ 
tion available from the seeding phase while annihilating the 
frequencies with indices k > K that are present in x but not 


Remark 2 In classical discrete-time signal processing, recov¬ 
ery of bandlimited signals is a two-step process. Firstly, a 
sparse regular signal whose non-zero values coincide with 
those of the signal to recover is generated. Secondly, the (zero) 
values not specihed in the sparse signal are extrapolated from 
the non-zero ones using a low-pass hlter. Our approach in 
this paper is to preserve this two-step methodology and use it 
to recover bandlimited graph signals. This provides a way to 
regenerate a desired signal in a graph - either estimated from 
samples or otherwise - by acting on a subset of (seeding) 
nodes. As it will be shown in the ensuing sections, for signals 
dehned on a general graph, the non-zero values of the sparse 
signal in the hrst step will not coincide with those of the 
signal to recover. This deviates from the classical concept of 
interpolation, which assumes that the non-zero values are the 
same than those of the signal to reconstruct. 

The practical advantage of studying recovery schemes that 
use graph hlters is twofold. First, they can be implemented 
distributedly, using only local exchanges among neighbors. 
Second, since graph filters can be used to model diffusion 
processes (e.g. the spread of an opinion in a social network), 
our results can be used to reconstruct signals in network 
applications that implement linear diffusion dynamics. 


III. Multiple node - single time seeding 

In multiple node - single time (MN-ST) seeding we consider 
the particular case where there is only t — 1 seeding signal s 
so that X = s [cf. 0]- Denoting by F the amount of nonzero 
values in s, we interpret MN-ST seeding as having P seeding 
nodes that inject a single value, while the remaining N — P 
nodes keep silent; see left and center panels in Figure [T] Dehne 
the signal injected by node i as Si and assume, without loss 
of generality, that the seeding nodes are the P hrst ones. We 
therefore dehne the P x 1 and x 1 seeding vectors as 

sp = [si,...,spf, ( 8 ) 

s = ...,0]^. (9) 

Then, given a bandlimited signal y = Vp-yp- [cf. Q], our 
goal is to design H and s such that 

y = Hs, (10) 


where H has the particular structure of a graph filter (cf. Sec 
tion 


II-Al. Exploiting the fact that y is bandlimited, it is 
reasonable to write ( [TOl i in the frequency domain. To do this, 
both sides of ( [T0| ) are left multiplied by V“^, which yields 

y = V-iHs = V-iVdiag(^h)V-is = diag('5'h)s, (11) 


where we used ([^ for the second equality. Utilizing the fact 
that the seeding signal s is sparse [cf. 0] we may write its 
frequency representation as 

s = V-^s = V-^EpSp, 


( 12 ) 
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where, we recall, Ep := [ei,ep] is a tall matrix collecting 
the P first canonical basis vectors of size iV x 1. By substi¬ 
tuting ( [T^ into •[ID, our goal of designing H and s such that 
y = Hs can be reformulated as designing h and sp such that 

y = diag(’®'h) V“^EpSp, (13) 

which is a bilinear system of N equations and L+P variables. 
Leveraging the sparsity of y [cf. Q], the system of N 
equations in ( [T3] l can be split into two 

fx = diag('5fh) V^^EpSp, (14) 

On-k = diag(’4'h) V”^EpSp, (15) 

where, we recall, Ep: := [ep-+i,epr] collects the last 
N — K canonical basis vectors and On-k denotes the 
{N — K) X 1 vector of all zeros. Note that the conditions in 
( [T5| | are the same for any itT-bandlimited signal. On the other 
hand, the conditions in ( [l4| ) depend on the specific signal to be 
interpolated. A natural approach is to use the filter coefficients 
h - which are related to the global behavior of the graph - to 
guarantee that ( |T5] l holds, while using the seeding signal sp 
to satisfy ( [l4| ) and, hence, to guarantee that the output of the 
interpolation is y. In this way, the filter coefficients h to be 
designed do not depend on the particular signal to reconstruct. 

The conditions under which the mentioned approach is 
guaranteed to find a feasible solution are given in the form 
of two propositions. Ensuing discussions describe the actual 
procedure to interpolate the signal. 

Proposition 1 IfL > D (cf. Section [70] ), there exist infinitely 
many nonzero L x 1 vectors h* such that, after setting h = h*, 
is satisfied for any V ^ and Sp. 

Proof: Since ( [T5| ) has to hold for any seeding signal sp, we 
need El^tPh = 0. This requires h to belong to the kernel of 
the {N — K) X L matrix E^\P. Since E^\P is a Vandermonde 
matrix, its number of linearly independent rows is equal to the 
number of distinct eigenvalues in {^k}k=K+i^ which is D. 
Thus, the existence of a solution h* 0 requires L > D. ■ 

For L > D, E]^’®' is rank deficient and the dimension of 
its kernel space is L — D. Hence, setting h to any nonzero 
element of the kernel space will satisfy ( [TS] ). In what follows, 
we will assume that L = D + 1 and set the coefficients h* 
that solve ( fTSj l as the unit vector spanning the unidimensional 
kernel space of E]^’®'. For the case where all the eigenvalues 
{^k}LK+l are distinct, this implies that L = N — K + 1. 

Once the coefficients of the filter are designed, the next step 
is to find the optimum seeding signal. With h.*j^ := E]^ ’®'h* 
denoting the frequency response of the low-pass filter in the 
active frequencies, substituting h = h* into ( [T4l l yields 

y^ = diagCh^f) E^V-iEpSp. (16) 

For which the following result holds. 

Proposition 2 The system of K equations in is guaran¬ 
teed to have a solution with respect to Sp if the following two 
conditions hold: 

i) Afci f Xk 2 for all (Aki, A/c^) such that ki < K and k 2 > K, 

ii) rank(E^V-iEp) > K. 


Proof : Condition i) is required to guarantee that all the 
elements of are nonzero. We prove this by contradiction. 
Recall that the following facts holds true (cf. Proposition [^i: 
a) E^^h* = Ox_x', b) h* f Op and c) rank of E^®' is 
L—1. Assume, without loss of generality, that the element of 
hjf that is zero is the if-th one. Then, we can use a) to write 
= Ox-K+i- Condition i) and fact c) guarantee that 
has rank L; then, satisfying E|^_j^’®'h* = Ox-k+i 
requires h* = Op, which contradicts b). Hence, all the 
elements of are nonzero. This guarantees that diag(h^) 
is invertible, so that can be written as 


diag(hJf)-iyK = (E^V-iEp)sp, (17) 

where E^V“^Ep isa KxP submatrix of V~^. To guarantee 
that the system of equations in ( [T7| ) has at least one solution, 
we need condition ii) to hold. ■ 


Different from the time domain, where all the eigenvalues 
of S (frequencies) are distinct, in the more general graph 
domain there can be graph topologies that give rise to S with 
repeated eigenvalues. Condition i) is required because a graph 
filter H always produces the same frequency response if the 
corresponding eigenvalues are the same. Therefore, it is not 
possible for H to eliminate one of the frequencies without 
eliminating the other. An alternative to bypass this problem 


is discussed in Section III-A Condition ii) requires the rank 
of the KxP matrix E]fV“^Ep being at least K. At the 
very least, this requires P, the number of seeding nodes, to be 
equal to K, the number of frequencies present in y. For the 
particular case of P = K, if the conditions in the proposition 
are satisfied, the sp that recovers y is 


Sp = (E^V-iEp)-idiag(hJf)-iyK. 


( 18 ) 


However, if condition ii) is not satisfied, there may be cases 
where setting P — K fails. To see why this is true, notice that 
[V-i \k,p can be viewed as how strongly node p expresses 
frequency k. Suppose for example that there exists a k such 
that p = 0 for all nodes p — 1,... ,P, then it is not 

possible to reconstruct a signal y with yp ^ 0 using that set 
of nodes. This problem is also present when sampling graph 
signals by observing the value of the signal in a subset of 
nodes Q. 

Proposition [^ states conditions for perfect reconstruction 
in a noiseless setting. In noisy scenarios, the specific set of 
nodes selected to inject the seeding signal has an impact on 


the reconstruction error. This is analyzed in Section VI 


A. Filter degree reduction in MN-ST seeding 

The MN-ST reconstruction scheme requires a low-pass filter 
H of degree D, which grows with the size of the graph. 
Since the degree of H corresponds to the number of local 
interactions needed to implement the hlter, the communication 
overhead can be a problem for large graphs. In this context, we 
look for solutions that reduce the degree of H by increasing 
the number of seeding nodes P. This can be done by splitting 
the system of equations in ( [T3] l as [cf. ([T4|l-([T5]l] 

[y^, Ol-xf = E? diag(®/h) V-^EpSp, (19) 
On-p = Ep diag('®'h) V~^EpSp. (20) 
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Fig. 1: Two different schemes to reconstruct signal y. On the left, MN-ST seeding injects P = 2 values simultaneously (blue arrows), after 
which the low-pass filter H is applied. On the right, SN-MT seeding first injects a single value that percolates to the two neighboring nodes. 
After a second value injection at the same node, filter H completes the reconstruction. 


The filter coefficients must be obtained now to annihilate the 
N — P frequencies in ( |20l l and the P seeding nodes must 
inject a signal whose spectrum, after being filtered by H, 
matches that of the desired signal ( [T^ . Notice that ([T^-(|20|) 
can also be used when d -([T5| | fail due to a violation of 
condition i) in Proposition |2 More specifically, for every 
frequency index k 2 > K with the same eigenvalue as a 
frequency index ki < K we can induce a zero frequency 
coefficient in the reconstructed signal via the seeding values 
[cf. ( [T^ ] instead of through the low-pass filter [cf. ( |20| )] and, 
hence, drop condition tj as a requirement for recovery. Further 
notice that for ( |20l i to hold for any sp, the degree of the filter 
needs to be at least equal to the number of distinct eigenvalues 
in {^k}^=p+i (cf. Proposition ll|. In the extreme case of 
P — N, the trivial solution h = ^, 0,..., 0]^ (0-order filter) 
and Sp = y satisfies (|T9l)-(|20|. 


B. Relation to classical interpolation 

In the classical time domain, sine (low-pass) interpolation 
of a bandlimited signal leads to perfect reconstruction. If 
the sampling is performed at the minimum possible rate, the 
bandwidth of the low-pass filter has to be exactly the same 
than that of y. By contrast, if the signal is oversampled, the 
bandwidth can be larger. Equivalently, if more samples than 
the minimum required number are available, then the low-pass 
filter does not have to cancel all the frequencies not present 


in y. The analysis in Section III-A reveals that this is also the 
case when signals are defined in more general graph domains. 

The main differences between the MN-ST reconstruction 
scheme and classical time interpolation come from the fact 
that the basis V of a general graph shift S is not as structured 
as the Fourier basis F. A difference of particular relevance is 
that, for general graphs, the seeding values sp do not coincide 
with the values of the desired signal y. This contrasts with 
the classical interpolation of uniformly sampled time-varying 
signals, where sp is a subset of the signal y. In fact, it 
can be rigorously shown that requiring such a condition for 
general graphs would lead to an infeasible interpolation. To 
be concrete, suppose that P = K, so that s = [cf. 

and that sp- is equal to the first entries of y. We can then 
leverage the fact that s and y are sparse to write 

SK = E^y = E^Vy = E^VE^y^. (21) 

Secondly, we write the goal of y = Hs into the frequency 
domain as y = diag(h)V~^s and use again the sparsity of s 
and y to write 

yi<r = EKdiag(h)V"^Ep:Sif =E^diag(h)Ep-E]^V^EifSp- 

= diag(h^)E^V-iEKSK, (22) 


where hp- := E^h contains the first K components of h. 
Substituting ( |2T| l into ( |22| ) yields 

Yk = diag(hK)E^V-iE^E^VE^fy^. (23) 

Since ( [23] l must hold for all yp-, it can only be satisfied 
if diag(hp')E|^V“^Ep-E]^VEp- = I. This requires matrix 
(E^V'^Ep-E^VEp:) to be diagonal. While this is true when 
K = N, it is not true for a general K. However, in the 
time domain where V = F (see Remark [TJ, for some cases 
the multiplication of submatrices of F is guaranteed to be 
diagonal. For example, if the K seeding nodes are chosen 
uniformly (equally) spaced, then (E^V“^Ep-E]^VEp-) = 
K/NI. This implies not only that ( |23| l is satisfied, but also 
that all the entries in hp- must be set to N/K. In other words, 
the optimal low-pass interpolator after uniform sampling in 
the time domain has the same response for all the active 
frequencies, as known from classical signal processing. 

IV. Single node - multiple time seeding 

In single node - multiple time (SN-MT) seeding, we con¬ 
sider the particular case where all the t = P seeding values 
are injected at a single node; see right and center panels in 
Figure [T] To be more specific, assume without loss of gener¬ 
ality that the first node is the one injecting the seeding values, 
so that the seeding signal at time t is of the form = 
0,..., 0]^. Then, define sp := ..., to be 

a P X 1 vector grouping the seeding values. We present the 
relation between the seeding values sp and the output of the 
seeding phase x in the following lemma. 

Lemma 1 The frequency representation of the intermediate 
signal X in SN-MT seeding is given by 

X = diag(ei)4'sp, (24) 

where ei := V"^ei is the frequency representation of the first 
canonical basis vector. 

Proof: Since x is obtained after P injections of seeding values 
following the dynamics in Q, it holds that 

X = x(^-i) = Er=”o' S's(^-i-0ei. 

(25) 

Equation ( |25] l relates the signal x to the successive inputs 
of the seeding node and can be interpreted as the application 
of the graph filter 


(26) 
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of degree P — 1 to the canonical basis vector ei. Building on 
this interpretation, we may use 0 to write 

X = diag(’5'sp)ei, (27) 


and, by exploiting the fact that for generic vectors a and b it 
holds that diag(a)b = diag(b)a, the lemma follows. ■ 


The proof of Lemma exploits the reinterpretation of the 
seeding phase as the application of a graph filter H, whose 
coefficients are given by the seeding values, to the canonical 
basis vector ei. Equation ( |24| l reveals how x depends on the 
structure of the graph and the seeding values sp, as well as 
on the particular node chosen to inject the values via ei, whose 
elements represent how strongly the node expresses each of 
the graph frequencies. 

The next step is to analyze the output of the filtering 
phase in the frequency domain z. To do this, recall that h* 
denotes the coefficients of a low-pass filter (cf. Section II-A[ ) 
that eliminates all frequencies with indices k > K. Defining 
h* := ’J'h*, we may analyze the application of the low-pass 
filter in the frequency domain as 


zk = E^diag(h*)x = E^diag(h*)EKE^x. 


(28) 


Further recalling that := E^h* and substituting ( |^ into 
( |28l l, we obtain [cf. ([H)] 

Yk = diag(h'^)E^diag(ei)’4'sp. (29) 


Expression ( [29l l relates the frequencies present in y to the 
seeding values sp. Provided that K < P, the following 
proposition states the conditions under which ( |29l l can be 
solved with respect to sp. 


Proposition 3 Let Ui be the number of values in {^i\k\k=i 
that are zero and let Di be the number of repeated values 
in Then, the system of K equations in ( |29| l is 

guaranteed to have a solution with respect to Sp if the 
following two conditions hold: 

i) Afcj f Xk 2 for all {\k ^, ) such that ki < K and k 2 > K, 

ii) Ui =0 and Di = 0. 

Proof; If we rewrite ( |29l l as 

Yk = (diag(h5f)) (E^diag(ei)Ep:) (E|^^) sp, (30) 

then it becomes clear that conditions i) and ii) ensure in- 
vertibility of the two square matrices, and full row rank of 
the rectangular matrix E|^'i'. To be specific, condition i) 
is required to guarantee that all the entries of vector 
are nonzero and, hence, matrix diag(hjf) is invertible (cf. 
proof of Proposition]^. Condition t/i = 0 in ii) ensures that 
E^diag(ei)Ep- is invertible since it is a diagonal matrix with 
no zero elements in its diagonal. Finally, Di = 0 guarantees 
that E^^f has rank K whenever K < P since it is a row-wise 
Vandermonde matrix with no repeated rows. ■ 

Condition i), which is equivalent to that in Proposition]^ 
guarantees that the low-pass filter with coefficients h* does not 
eliminate any of the frequencies present in y. Condition ii) 
states requirements for recovery on both the seeding node 
and the global structure of the graph. The seeding node is 
required to be able to act on every active frequency (Ui = 0), 
while the graph is required to have every active frequency 


distinguishable from each other (Di = 0). Condition ii) 
ensures that the rank of matrix E]^diag(ei)’®' is equal to 
K when P > K, guaranteeing that ( ]29| l can be solved with 
respect to sp. For the particular case of P = K the seeding 
values can be found as 

Sp = (E^diag(ei)'a') ^ diag(h;^)“^yK. (31) 

When comparing the conditions ii) in Propositions ]^ and ]^ 
we observe that for MN-ST seeding we should require a rank 
condition on a submatrix of V“^. By contrast, for SN-MT 
seeding, the Vandermonde structure of \P allows reformulating 
the rank condition in terms of the graph related quantities Ui 
and Di, providing further insight on specifying the situations 
when recovery is possible. This dual behavior is also present 
when sampling graph signals. When following “selection 
sampling” |]^, which is the counterpart of MN-ST interpo¬ 
lation, perfect reconstruction depends on the invertibility of 
a submatrix of V, whereas when following an “aggregation 
sampling” scheme Q, which is the counterpart of SN-MT 
interpolation, the conditions for perfect reconstruction can be 
written in terms of specific graph related quantities. 

Even though Proposition ]^ guarantees perfect recovery 
under SN-MT seeding in a noiseless case, in noisy scenarios 
the selection of the seeding node is essential to reduce the 
reconstruction error. This is analyzed in Section ]Vl] under a 
more general seeding scheme. 


A. Filter degree reduction in SN-MT seeding 

Mimicking the filter degree reduction technique presented 
in Section ]TII-A| SN-MT seeding can also implement a lower- 
degree filter if a higher number of seeding values is injected. 
To achieve this, we need the additional seeding values to 
generate a signal whose spectrum is zero for the inactive fre¬ 
quencies that are not eliminated by the filter. More specifically, 
the seeding values sp and the filter coefficients h have to 
satisfy [cf 

[y^, = diag(Ep’S'h)Epdiag(ei)’®'sp, (32) 

Ojv_p = diag(Ep’®'h)Epdiag(ei)’5'Sp, (33) 


where N — P is the number of frequency coefficients elim¬ 
inated by the low-pass filter h. As done in Section III-A h 


will be designed to solve ( ]3^ for any choice of sp, while sp 
will be chosen to solve the P equations in ( ]32| l. A sufficient 
degree for h is presented next. 


Proposition 4 Let U 2 be the number of values 
in that are zero and D 2 be the 

number of repeated values in {Xk\k^Ku’ where 
ICu := {k \ K < k < N and [eij^ ^ 0}. Then, ( ]3^ 
can be solved with respect to h/or any choice of sp provided 
that L—1 > max{0, N — P—U 2 — 02 ). 

Proof: Notice that in ( ]3^ the filter eliminates the last N — 
P frequencies, however, since the ordering is arbitrary, any 
subset of N — P frequencies (not containing the K first ones) 
can be chosen to be annihilated by h. Thus, our objective 
it to show that the proposed degree is enough to nullify a 
particular choice of N — P frequency coefficients. Define as 
TZ the set of indices corresponding to zero elements in ei or 
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repeated rows in Since condition ii) in Proposition must 
be satisfied - otherwise, perfect recovery would be infeasible 

the cardinality of TZ is U 2 + D 2 and every index in TZ must 
be greater than K. 

First assume that U 2 + D 2 < N — P and pick the 
N — P frequencies to be eliminated to include the ones in 
TZ. This is equivalent to picking a frequency ordering such 
that every index in TZ is greater than P. Based on TZ, define 
the selection matrices := ..., for all 

ki € TZ and E-ji := [efe^, ..., for all 

ki G {P +1,..., N} \ TZ where \ represents the set difference 
operator. Hence, the system of equations in ( [3^ can be split 
into two 

ON-P-U 2 -D 2 = diag(E^^h)E^diag(ei)«'sp, (34) 

OU 2 +D 2 = diag(E^^h)E^diag(ei)«'sp. (35) 

Condition ( |34l i can be guaranteed for any sp if h = h*, where 
h* are the coefficients of a low-pass filter of degree L—1 = 
N — P —U 2 — D 2 , as stated by the proposition. To complete 
the proof, we need to show that h = h* also guarantees that 
© holds. To see why this is the case, notice that U 2 rows 
of Ei^diag(ei)’®' are exactly zero, trivially satisfying ( |T5] l for 
any sp. Also, each of the remaining D 2 equations in ( |T5| ) 
corresponds to a repeated eigenvalue and, thus, can be obtained 
by multiplying one of the N — K — U 2 — D 2 homogenous 
equations in ( |34l i and ( [32l i by a scalar, guaranteeing that h* 
also solves these D 2 equations. 

For the case where U 2 + D 2 > N — P, we pick a frequency 
ordering such that every index greater than P is contained in 
TZ. Thus, ( [33| is implied by the homogenous equations in ( [32l i, 
and no filter (degree 0) is needed. ■ 

Proposition explicitly states that every additional seeding 
value decreases the required filter degree. However, in contrast 
to the situation for MN-ST, this reduction of the filter degree 
does not entail a reduction in the number of applications 
of the graph-shift operator, because it requires the length 
of the seeding phase to be extended. More interestingly, the 
additional seeding values can be used as a mean to guarantee 
perfect reconstruction when condition i) in Proposition[^is not 
satisfied, as explained in Section III-A| for MN-ST seeding. 

When P > N—U 2 —D 2 the seeding phase suffices to recover 
the signal. This can be of interest in scenarios where the graph- 
shift operator S describes an intrinsic graph diffusion dynamic 
and the design of the filter coefficients is not feasible. It is also 
of interest if y is not bandlimited. See Section m for further 
discussions. 

B. Relation to classical interpolation 

When S = A^c, applying the Zth power of S to a signal 
amounts to shifting the signal I time instants. Consequently, 
the intermediate signal x obtained after the seeding phase in 
SN-MT reconstruction coincides with the seeding signal s in 
MN-ST reconstruction, provided that the seeding nodes are 
chosen adjacent to each other. Moreover, for the extreme case 
of the number of seeding values P being enough to eliminate 
the filtering phase, which entails L — 1 = 0, Proposition 
requires setting P = N, because both U 2 and D 2 are zero if 
S = Aiic. The design of the P = N seeding values Sp that 
guarantee that x = y can be carried out trivially by setting 
Sp = [a:i,... ,a;Ar] = y. 


V. Multiple node - multiple time seeding 

In the more general multiple node - multiple time (MN-MT) 
seeding scheme, we can have several seeding signals (t > 1) 
and we do not assume any structure on so that any node 
may inject a seeding value at any given time. We concatenate 
the r seeding signals into the Nt x 1 vector s defined as 
s :=vec([s("-i),s("-2),...,sW]^). 

Defining the NxN'^ matrix © := [diag(ei),... ,diag(ejv)], 
we may relate x to s as stated in the following lemma. 


Lemma 2 The frequency representation of the intermediate 
signal X in MN-MT seeding is given by 

X = 0(10 ’S')s, (36) 


where ® represents the Kronecker product. 


Proof: If we denote by x := x*^"^ the signal obtained after 
the seeding phase, it holds that [cf 0] 




1=0 


1=0 



N 


= Eh. 


2=1 


where the filter is given by 


r-l 


H. = E 


cl 


(37) 


(38) 


1=0 


Writing the input-output relationship of those filters in the 
frequency domain, we have that [cf. 0] 


N 


2=1 


N 

X = y^diag(’®'Sj)ei = diag(e^)^s. 

2 = 1 


(39) 


Recalling the definitions of 0 and s, the sum in ( [39| ) can be 
written in matrix form, giving rise to ( [36] l. ■ 

As was the case for Lemma[^in SN-MT seeding, Lemma]^ 
leverages the reinterpretation of the seeding phase as the 
application of a filter - in this case, N different filters, one per 
node - to the canonical basis vectors [cf. ([37|]. The coefficients 
of the filter associated with the z-th node are given by the 
values injected by that z-th node [cf. (|38]l]. Notice that, as 
expected, ( |3^ reduces to whenever the seeding values 
are forced to be zero for every seeding node except for the 
first one. 

To analyze the output of the filtering phase z, recall that 
hjf = E](^’®'h* denotes the response of a low-pass filter in the 


active frequencies. Mimicking the procedure in Section IV 


we 


find that the active frequency coefficients in y can be written 
in terms of the seeding values s as 


yji = zk = diag(h;^)E^©(I (g) 'i'js. 


(40) 


The system of equations in ( |40l i is underdetermined, since the 
K values in yp- can be reconstructed using the Nt values in 
s. However, our focus is on the case where only P Nt 
seeding values are injected during the seeding phase. To this 
extent, we introduce the P x Nt selection matrix C whose 
elements are binary Cij G {0,1} and satisfy JJj Cij = 1 and 
Ei C'ij < 1 for all z and j, respectively. Since the matrix has 
exactly one 1 in every row, C selects P seeding values among 
the Nt node-time pairs. If we denote by Sp := Cs the vector 






containing these P seeding values, (|4gi can be rewritten as 
[cf. @ and @] 

Yk = diag(h;,)E^©(I 0 ’5')C^Sp. (41) 

To resemble the structure of previous sections, the conditions 
under which ( |4T] l can be solved with respect to Sp are given 
in the form of a proposition. 


Proposition 5 The system of K equations in ED is guaran¬ 
teed to have a solution with respect to Sp if the following two 
conditions hold: 

i) ^ki 7 ^ A/c 2 for all (Afc^, \k 2 ) ^uch that ki < K and k 2 > K, 

ii) rank(E^©(l0 > K. 

Proof: Condition i) is required to guarantee that diag(hj^) is 
invertible (cf. proof of Proposition]^. This allows us to rewrite 
( |4T] l as 

diag(h'Jf = E^©(I 0 '5f)C^Sp. (42) 

To guarantee that the system of equations in ( |4^ has at least 
one solution, we need condition ii) to hold. ■ 

Condition i), also present in Propositions]^ andguarantees 
that the filtering phase does not annihilate any of the frequen¬ 
cies present in y. Condition ii) requires, at the very least, 
P > K. However, there may be cases where setting P = K 
can fail as stated in the discussion ensuing Proposition ^ 
Mimicking the developments in Sections |III-A| and IV-A| 
for the general case of MN-MT seeding, additional seeding 
values can be used to reduce the degree of the low-pass filter. 
Indeed, for every extra seeding value the degree of the filter 
needed decreases by one, reducing the communication cost 
of the reconstruction scheme. Moreover, these extra seeding 
values can be used to obtain perfect reconstruction even when 
condition i) in Proposition ]^ is violated, as explained in 
Section IIII-AI 

The selection matrix C can be designed so that condition ii) 
in Proposition ]^ is satisfied, guaranteeing perfect recovery. 
Furthermore, for the cases in which perfect reconstruction is 
infeasible due to, e.g., the presence of noise, the choice of C 
can be optimized to achieve robust recovery, as analyzed in 
the following section. 


Remark 3 When C in ( ]4T] i selects the first P elements of 
s, iD reduces to ( ]29l l and SN-MT reconstruction is recov¬ 
ered. Similarly, if C selects the elements of s in positions 
1, r -f 1,..., Pt -f 1, then ( ]4T] i reduces to ( fTh] ) as in MN-ST 
reconstruction. 


VI. Imperfect reconstruction 


We study two settings where perfect reconstruction is infea¬ 
sible; insufficient number of seeding values (Section VI-A| l and 
additive noise in the injections (Section VI-B i. The analysis 
is focused on the MN-MT seeding scheme, since the results 
for MN-ST and SN-MT can be obtained by particularizing the 
value of the selection matrix C (cf Remark]^. 


A. Insufficient seeding values 

When the number of seeding values P is not enough to 
achieve perfect reconstruction, the goal is to minimize a pre¬ 
specified error metric between the reconstructed signal z and 
the original AT-bandlimited graph signal y. Three different 
design scenarios are considered. In the first one, the seeding 
values Sp are designed assuming that both h and C are fixed. 
The second scenario addresses the joint design of Sp and h. 
In the last one, the joint design of Sp and C is performed. 

1) Designing the seeding values Sp: Assume that condi¬ 
tion i) in Proposition ]^ holds and recall that h* stands for 
the coefficients of a low-pass filter that eliminates all the 
frequencies k > K. Then, the first K frequency coefficients 
Zp; of the reconstructed signal z are obtained as [cf. ED] 

Zp = diag(hjf )E]^©(I 0 ’®')C^Sp. (43) 

Since we assume insufficient seeding values, i.e., P < K, 
obtaining yp = zp is in general infeasible. A reasonable 
approach is to design Sp to minimize the energy of the 
reconstruction error. Defining the matrix 

:= diag(hJf)E^©(I 0 ^), (44) 

the optimal seeding values Sp can be obtained as 

Sp := argmin||y-Vif$ifC^Sp|| 2 , (45) 

§.p 

where, we recall, Yp := VEp-. The minimization problem 
in ( ]45| ) has the well-known closed-form solution | ]24l 

s*p = Vf y, (46) 

where we assume that the fixed seeding locations C lead to 
a matrix ^pC'^ that has full column rank. With e y — z 
denoting the reconstruction error, its energy can be written as 

l|e |!2 = y"V;f (l - Vf y. 

(47) 

Notice that, since h* is given, the reconstruction error is zero 
for the frequency components k > K. 

2) Designing the seeding values Sp and the filter coeffi¬ 
cients h.- When perfect reconstruction is infeasible, carrying 
out a separate optimization of h and Sp, where h is designed 
to filter the frequencies not present in y and Sp is designed 
to match the spectrum of y in the active frequencies, is not 
jointly optimal. Minimization of the reconstruction error by 
jointly designing Sp and h is briefly discussed next. Notice 
that the N frequency coefficients - as opposed to just the first 
K coefficients - of the reconstructed signal z are [cf. ED] 

z = diag(’®'h)E^©(l0'5')C^Sp. (48) 

Hence, if the objective is to minimize llejlf, we have that 
{sp,h*} := argmin||y - Vz ||2 (49) 

{s^,h} 

= argmin ||y - Vdiag('5'h)E|^©(I 0 Sp\\l, 

{sp.h} 

which is a bilinear optimization. Bilinear problems are non- 
convex, but there is a large amount of works dealing with their 
analysis and efficient solution l@-|]23. 

The formulation in ( ]49l l considers that h* can be chosen as a 
function of the signal to reconstruct y. In applications where 
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this is not feasible, formulating the optimal design requires 
additional assumptions on y. If the distribution of y is known, 
a two-stage stochastic programming approach can be pursued 
In the second stage, h in ( |49| ) is considered given and 
the optimal Sp is obtained as the minimizer of ||e(h, y, Sp)|||, 
which is a function of h and y [cf. (|46|)]. In the first stage, the 
solution of the second stage Sp(h,y) and the distribution of 
y are leveraged to write the expectation of the reconstruction 
eiTor in ( |49l l as e(h) Ey[||e(h,y,Sp(h,y))|||], which only 
depends on h. The optimum h* is then the minimizer of the 
expected error e(h). Notice that this two-stage approach is 
used in Sections and |V] to find conditions for perfect 

recovery, where bandlimitedness is the prior knowledge of y. 

3) Designing the seeding values Sp and the seeding loca¬ 
tions C; Suppose now that one can select the specific nodes 
and time instants where the injections take place. This amounts 
to choosing the P entries of C that are non-zero, which is a 
combinatorial problem. Although for small networks one could 
try all possible choices of C and select the one leading to 
the smallest reconstruction error, for general networks a more 
scalable approach is required. To formulate the corresponding 
optimization problem Sp = C s is substituted into ( |45] l. After 
that, the product C^C is rewritten as diag(c) where c is a 
binary selection vector of dimension Ntx\. Note that having 
Ci = 1 indicates that at time t = modr{NT — i) the node 
{i -f f)/r injects a seeding value. With this notation, the joint 
design of Sp and c amounts to solving 

{s*,c*} := argmin||y - V/f$xdiag(c) s|if -f 7 ||c||o 

{s,c} 

s.t. ce{0,l}^^, (50) 


where is defined in ( |44l l. In ( |50l l each seeding location 
used is penalized with a constant cost 7 . By tuning 7 , the 
desired level of sparsity of c can be achieved. Problem ( |50l l can 
be further simplified by setting d := diag(c)s and requiring 
sparsity on d 

d* := argmin ||y - Yk^k^^WI + 7 ||d||o. (51) 

d 


Among other advantages, the formulation in ( [5T| l is amenable 
to relaxations that reduce the computational complexity re¬ 
quired to find a solution. A straightforward approach is to 
relax the problem by replacing the 0 -norm with the 1 -norm to 
obtain a convex formulation. 

As in problem j4^, the design in ( |50l l and its subsequent 
simplification in ( |51| l assume that the seeding nodes can be 
chosen as a function of y. For applications where this is 
not convenient, a two-stage stochastic programming approach 
similar to the one described for ( |49] l can also be used in solving 
®. Last but not least, although computationally challenging, 
a joint optimization of Sp, h and C can be pursued by 
combining the approaches in Sections VI-A2 and VI-A3| 


applications of S. The goal of this section is to quantify the 
reconstruction error and to discuss seeding selection schemes 
tailored to these operating conditions. Their performance will 
be illustrated via numerical simulations in Section I VII I 
Let us assume that the injected signal is Sp-fwp, where wp 
is a P X 1 noise vector with zero mean and covariance R^. 
The active frequencies of the reconstructed signal can then be 
written as 2.x = -f Wp) [cf. ( |4T] i and (|44li]. From 

this, we may obtain the reconstruction error as 

e = Vif (zif — yif) = Vx^icC^wp, (52) 
with covariance matrix 

R, = E(ee^) = C^RwC^f Vf. (53) 

Ideally, C should be designed to select the seeding nodes and 
time instants that minimize the reconstruction error, which 
can be quantified as a function of R^. In what follows, we 
will focus on minimizing the mean squared error (MSE), 
which is achieved by minimizing trace(Re). However, similar 
approaches can be followed to minimize other commonly used 
error metrics such as Aniax(R-e) and log(det(Re)) p^ . 

To illustrate the design of C, two particular scenarios are 
considered. In the first one, we assume i) that the frequency 
coefficients yx are zero mean with covariance Ry = I and 
ii) that the noise wp is zero mean with covariance R^ = 
Edlsplll)!. Note that assumption ii) is meaningful if the 
system operates under a constant signal-to-noise ratio (SNR) 
regime. In the second scenario, we assume the noise is also 
uncorrelated but its power is independent of that of the seeding 
signal, so that Rw = ct^I. 

For convenience, the optimal seeding strategy for the first 
scenario is presented in the form of a lemma. 


Lemma 3 Suppose that yx and wp are drawn from zero- 
mean distributions with covariances Ry = I and R^ = 
Edisplll)!, respectively. Then, the selection c* that mini¬ 
mizes the MSE of the reconstruction is given by 

c*:=argmin trace((4>p-diag(c)$^) ^)trace('4>p'diag(c)$^^ 

s.t. cG{0,17", ||c||o = P. (54) 

Proof : To prove the lemma, we need to show that the 
minimization of the objective in ( |54l i is equivalent to the min¬ 
imization of trace(RE). By substituting R^ = cr^ Ed|sp|j|)I 
and diag(c) := C^C into ( |5^ , it follows that 

trace(Re) = Ed|sp|| 2 )trace(Vp-$p'diag(c) 4 >^Vj^). 

(55) 

Since the trace is invariant to cyclic permutations and 
V^Vp- = I, we have that 


B. Noise when injecting the seeding values 

The conditions for perfect reconstruction stated in Propo¬ 
sitions mil] and [5| re quire the seeding values to be the exact 
solution of ( |16| l, ^9] ), and iD’ respectively. However, in real 
applications, the injected values can be corrupted with additive 
noise. This noise can be either attenuated or amplified when 
the signal percolates through the graph via the successive 


trace(Re) = cr^ Ed|sp||2)trace(4>p:diag(c)$2). (56) 

To find an expression for Ed|sp|||), we leverage the fact 
that yx = ^ifC^Sp [cf. ED and ED] to write Sp = 
C$ 2 ($p-C^C$^)“^yp- and, consequently, to write ||sp ||2 
as 

|]sp||i = sfsp = yf ($p-diag(c)$f )-iyK. 


(57) 






10 





1 

0.31 

-0.44 

-0.15 

-0.22 

-0.04 

0.12 

0.54 

-0.27 


1 

-0.03 

-0.00 

-0.00 

-0.00 -0.01 

-0.03 

-0.08 

-0.18 

2 

-0.42 

0.86 

0.44 0.52 

0.49 

0.22 

0.02 

-0.49 


2 

-0.02 

0.01 

-0.01 

0.01 

-0.02 

0.04 

-0.10 

0.31 

3 


0 

-0.43 -0.53 

-0.65 

-0.75 

-0.73 

-0.38 

OJ 

3 

-0.03 

0.07 

-0.06 

0.06 -0.09 

0.18 

-0.44 

1.24 

X A 
0) ^ 

0 

0.31 

-0.01 

-0.07 

-0.14 

0.02 

0.21 

0.39 

c 

4 

0.16 

0.05 

0.15 

0.40 0.82 

1.27 

1.36 

1.03 

c 5 

0 

-0.42 

0.07 0.17 

0.27 

-0.07 

-0.05 

-0.82 


5 

0.29 


0 

0 

0 

0 

0 

0 

<1) 6 

0 

0 

0.08 0.09 

0.04 

0.16 

0.16 

0.37 

c 

0) 

6 

-0.07 

0.19 

0.02 

0 

0 

0 

0 

0 

0 7 


-0.11 

0.14 0.42 

0.59 

0.75 

0.30 

0.76 

S' 

7 

0.25 

-0.22 

-0.45 

-0.85 

-1.11 

•0.87 

-0.28 

0 

^ 8 


0.31 

-0.10 -0.04 

0.22 

0.47 

0.57 

0.31 

£ 

8 

0.18 

-0.11 

-0.18 

-0.29 -0.29 

-0.13 

0 

0 

g 


0 

-0.09 -0.46 

-0.75 

-0.85 

-0.40 

-0.77 


9 

-0.21 

0.53 

0.37 

0.22 

0 

0 

0 

0 

10 

0 

-0.11 

-0.11 

-0.24 

-0.53 

-0.49 

-0.81 

0.22 


10 

0.11 

-0.20 

-0.24 

-0.26 -0.13 

0 

0 

0 


0 


2 

3 

4 

5 

6 

7 



0 


2 

3 

4 

5 

6 

7 


Time 


Time 


(a) (b) (c) 

Fig. 2: Perfect recovery of a bandlimited graph signal, (a) The graph Q, the target signal to recover y and its frequency representation y. 
(b) Evolution of the reconstructed signal. The seeding and filtering phases are separated by a dotted line and the recovered signal is framed 
in red. (c) Evolution of the frequency components of the reconstructed signal. Successive annihilation during the filtering phase is observed. 


Using the expression for the expected value of a quadratic 
form, it follows that 

E(l|sp|l 2 ) = trace(^($ifdiag(c)$^)"^^ . (58) 

Upon replacing ( |58| ) into ( |5^ and recalling that does not 
depend on c, the expression in ( |54l l follows. ■ 

The statistical assumption on yx allows us to design c* such 
that the expected performance of the reconstruction scheme 
is optimized. In this way, the choice of the seeding nodes 
and instants is independent of the particular signal being 
reconstructed. 

Even though obtaining general relaxations to efficiently 
approximate the non-convex problem in ( |54| ) is out of the 
scope of the paper, we can gain intuition by specializing 0 
for time-varying signals, i.e., by setting S = A^c- For SN- 
MT seeding, where designing c boils down to selecting the 
seeding node, it can be shown that the objective in 0 does 
not depend on the particular node chosen. This is as it should 
be, since in the directed cycle every node is topologically 
indistinguishable from the others. For MN-ST seeding, the 
best strategy is to uniformly distribute the seeding nodes, as 
we formally state next. 

Proposition 6 Suppose that the problem in ( |54| i is particu¬ 
larized for the case of MN-ST seeding of time-varying signals 
using an ideal low-pass filter. Then, if K = P = N/9, it holds 
that the optimal seeding strategy selects the nodes in positions 

l,l + e,...,l + {K-l)9. 



MN-ST 

SN-MT 

MN-MT 

% of recovery 

91.8 

96.4 

94.4 

Min error 

.001 

.032 

.003 

Median error 

.048 

.349 

.066 


TABLE I: Recovery performance for the three seeding schemes. We 
restrict MN-MT to consist of two seeding nodes injecting two values. 

which does not depend on c. Hence, the optimal c* in ( |59l ) 
can be found as the one minimizing trace 

If we denote by the K eigenvalues of M, our goal 

is then to find the c* that minimizes ^Hi- Given that all 7 ^ 
are nonnegative (M is positive semi-definite) and ( [ 6 OI 1 implies 
that = K'^/N, the minimization is achieved by setting 

7 i = 72 = ... = 7 if = K/N. Hence, if we show that uniform 
seeding leads to 7 ^ = K/N for all i, the proof concludes. To 
show this, notice that under uniform sampling 

diag(c)FE;f = (61) 

where F^^^ is the Fourier basis of size K x K. Hence, M = 
K/Nl [cf ( |59| )] and every eigenvalue of M equals K/N. ■ 

In words, even though for the noiseless case any seeding 
selection strategy satisfying the conditions in Proposition is 
equally optimal, uniform seeding in directed cycles is the best 
MN-ST scheme when noise is present in Sp. 

A second scenario of interest are setups where the additive 
noise at different value injections is uncorrelated and of fixed 
power, i.e., = cr^I. In this case, 0 can be rewritten as 

R, = a2V^f$Kdiag(c)$f Vf. (62) 


Proof: When S = Adc we have that: a) V = F; b) = oIk 
for some constant a where Ip- is the K x 1 vector of all 
ones - since we are considering an ideal low-pass filter 
and c) c only can take nonzero values in positions i = l,l-|- 
T, ..., 1-|-(A—l)r (since we are considering MN-ST seeding). 
Leveraging a), b) and c), problem ( |5^ can be reformulated as 

c* := argmin trace (M~^) trace (M) (59) 

C 

s.t. M = E^F"diag(c)FEK, c e {0,1}"^, ||c||o = A, 


where c selects K seeding nodes out of the N possible ones. 
First, notice that trace(M) does not depend on the particular 
choice of c. To see why this is true, we denote by 1(c) the 
set containing the indices of the K seeding nodes selected by 
c. Then, we can exploit the structure in F to write 


K-l 

trace(M) = H 

iGX(c) j=0 


Vn 


K^ 
A ’ 


The design of c that minimizes the MSE of the reconstruction 
is the solution of the following linear integer program 

c* := argmin trace(R£) = argmin trace(€>p-diag(c)$)^) 

C C 

s.t. l|c||o = -P (63) 

which can be approximated by relaxing the binary and 0 -norm 
constraints. 

It turns out that the solution of ( | 6 ^ promotes the injection 
of seeding values at nodes that weakly express the active 
frequencies, i.e., nodes j such that the values for k < K 
are small. This occurs because the noise power is fixed and 
those nodes require the injection of seeding signals with high 
power, leading to a high SNR. 


VH. Numerical experiments 
We illustrate the reconstruction schemes in noiseless and 
noisy scenarios using synthetic (Section [VII-A[ ) and real-world 
graphs (Sections VII-B and VII-C| l. 













































II 


5 0-6 
-io.5 


© MNST 
•*-SNMT 
-•/i-MNMT 


2.5 3 3.5 

Seeding Values 


the reconstructed signal is sparse during the seeding phase, 
consisting of the hrst two time instants. More specihcally, for 
t = 0 the signal attains nonzero values only for the seeding 
nodes [cf. 0] and for t = 1 the signal remains zero for 
every node outside of the one-hop neighborhood of the seeding 
nodes. During the filtering phase - times t = 2 to t = 7- 
signal values are successively exchanged between neighboring 
nodes in order to hnally recover y at time t = 7. Figure 2c 
helps to understand the operation of the hltering phase. The 
signal X obtained after the seeding phase (t = 1) contains 
every frequency not active in the desired signal y. Thus, 
in every successive time instant, one of these frequencies is 
annihilated. E.g., at time t = 2 the frequency with index i = 5 
is eliminated and at f = 3 the frequency i = 6 is eliminated. 
In this way, at time t = 7 every frequency not active in y has 
been annihilated and perfect recovery is achieved. 

To compare the reconstruction performance of MN-ST, SN- 
MT, and MN-MT seeding, we generate 1000 Erdos-Renyi 
graphs with 10 nodes and edge probabilities between 0.2 
and 0.4. On each graph we dehne a 4-bandlimited signal 
and try to recover it through the three seeding schemes 
presented; see Table We restrict the MN-MT schemes to 
those consisting of two seeding nodes injecting two values 
each. We hrst compute the recovery percentage of the three 
schemes in noiseless scenarios. More specihcally, for a given 



Fig. 3: Reconstruction errors when recovering a signal in a social 
network with insufficient seeding values. 


A. Synthetic graph signals 


Eigure 2a represents a graph Q with iV = 10 nodes and 
adjacency matrix A generated using an Erdos-Renyi (ER) 
model with edge probability 0.3 pO) . Dehne the graph-shift 
operator S = A and let y be a signal to be recovered. Though 
seemingly random in the node domain, the structure of y is 
highly determined by Q. Indeed, y has bandwidth K = 4, 
as can be observed from its frequency representation y in 
Eigure 

The hrst set of experiments illustrates the perfect recovery 
of y when P = 4 seeding values are injected into Q followed 
by a low-pass hlter of degree N—P = 6. The reconstruction is 
carried out using MN-MT seeding (Section |V]l where nodes 1 
and 2 act as seeding nodes and each of them injects a seeding 
value for time instants t G {0,1}. After the seeding phase, a 
hlter that successively annihilates the N — K = Q frequencies 
not active in y is implemented [cf. 0]. The evolutions 
of the reconstructed signal and its frequency representation 
are depicted in Eigures and respectively. Notice that 
perfect reconstruction is achieved since the last column in 
both hgures coincide with y and y. Eigure 2b illustrates that 


Brain Regions 

Fig. 4: Heat map of the adjacency matrix A of brain graph Q. 

graph and signal to recover, we test for perfect recovery for 
every possible combination of seeding nodes. For example, 
there are 210 ways (10 choose 4) of selecting the seeding 
nodes in MN-ST while there are only 10 ways of selecting 
the single seeding node in SN-MT. If, e.g., 9 out of these 10 
ways lead to perfect recovery, then the recovery percentage 
for SN-MT on that particular graph is 90%. The values in 
Table [I] correspond to the averages of these percentages across 
the 1000 graphs generated. Notice that the highest recovery 
percentage of SN-MT suggests that condition ii) in Proposi¬ 
tion is more commonly satisfied in random ER graphs than 
the respective conditions in Propositions and We then 
introduce noise i n the injections following the constant SNR 
model in Section VI-B for cr = 10“^. Denoting by z the signal 


obtained from the reconstruction and by y the desired signal, 
we define the reconstruction error as e = |lz —y|j 2 /||y|| 2 - For 
every given graph and signal y, we record the minimum and 
median e for every possible choice of seeding nodes within 
each reconstruction scheme. In Table |I] we report the median 
of these values across the 1000 graphs generated. As it turns 
out, MN-ST is an order of magnitude more robust than SN-MT 
both in terms of minimum and median error. Finally, observe 
that MN-MT seeding presents an intermediate behavior both 
in terms of recovery percentage and reconstruction error. 

B. Influencing opinions in social networks 

Consider the well-known social network of Zachary’s karate 
club HD represented by a graph Q consisting of 34 nodes or 
members of the club and 78 undirected edges symbolizing 
friendships among members. Denoting by L the Laplacian 
of Q, define the graph shift operator S = I — aL with 
a = l/Aniax(L). A signal y on can be interpreted as 
a unidimensional opinion of each club member regarding 
a specihc topic, and each successive application of S can 
be seen as an opinion update influenced by neighboring 
individuals. Bandlimitedness of y implies that the opinion 
discrepancies between neighbors are small. In this context, 
signal reconstruction can be interpreted as the problem of 
inducing a desired global opinion prohle by influencing the 
opinion of a subset of members. 

We analyze the recovery performance when the number 
of seeding values is insufficient (Section |VI-A i. For this, we 
generate a signal y of bandwidth K = 5 and try to reconstruct 
it using P seeding values for P = 1,..., 5; see Figure 
For every P, we hnd the combination of seeding values and 
locations that minimizes the reconstruction error within each 
seeding scheme (cf. Section VI-A3|l. E.g., if P = 2 and we 
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Fig. 6: Inducing a brain state in the presence of noise, (a) Histogram of the reconstruction error for different choices of the six seeding 
nodes, (b) Frequency of appearance of each brain region among the configurations achieving the lowest reconstruction errors, (c) Anatomical 
location of the six regions most used in robust seeding configurations. 



Fig. 5: Initial yi (red) and target yt (blue) brain states. High activity is 
represented by positive activity levels while negative values represent 
low levels of activity. 


are analyzing SN-MT seeding, we consider every individual 
as possible seeding node and then choose the one achieving 
the minimum error. In Figure]^ we report the average of these 
minimum errors across 100 bandlimited signals. As expected, 
for P = 1 the three schemes coincide and for P = 5 perfect 
recovery is achieved for all of them. However, for intermediate 
values of P, MN-ST presents a considerably lower error than 
SN-MT. For P = 3 this implies that, when trying to induce 
a global opinion profile, it is more effective to influence the 
opinion of three individuals once than to influence the opinion 
of the same individual three times. The fact that MN-MT 
seeding presents the lowest reconstruction errors is expected 
since this scheme includes the other two as particular cases. 
Notice that in this case, as opposed to Table |I] we do not 
restrict MN-MT to the cases where multiple seeding nodes 
and multiple seeding values in each node are used. 

For the above analysis to hold true, we must be able to 
apply a low-pass Alter on the social network as required by the 
filtering phases of MN-ST, SN-MT, and MN-MT. This can be 
achieved by assuming that we can modify the rate of exchange 
of opinions in the network represented by a. Indeed, consider 
that after the seeding phase, the signal still percolates over the 
graph - people still communicate their opinions to neighbors 
- but we can modify the diffusion rate a; at each discrete time 
instant. Thus, after L — 1 interactions we obtain that 

z = nf=i^(I - a;L)x, (64) 

which is equivalent to applying an annihilating Alter [cf. 
to X. Notice that the Alter in ( |64| i is a polynomial on L rather 
than S. However, the frequency annihilation procedure is still 
valid since the eigenvectors - frequency basis - of L and S 
are equal. 


C. Inducing a brain state 


Upon dividing the human brain into the 66 regions of 
interest (ROIs) deflned in p2) , we build a weighted undirected 
graph Q whose nodes are the ROIs and whose edge weights 
are given by the density of anatomical connections between 
regions; see Figure]^ The flrst 33 ROIs are located on the right 
hemisphere of the brain while regions 34 to 66 correspond 
to their left counterparts. From Figure we see that most 
connections occur within the same cortical hemisphere with 
few inter hemispheric connections. We deflne the graph-shift 
operator S = A where A is the adjacency matrix of Q. 
The level of activity of each ROI can be represented by a 
graph signal y where larger values represent higher levels of 
activity. Successive applications of S on y model a linear 


evolution of the brain activity pattern |33|. As a method to 


inject seeding values to Q, we consider transcranial magnetic 
stimulation (TMS) | [T5| , a noninvasive method to stimulate 
ROIs. In this context, reconstructing a brain signal amounts 
to inducing a speciflc brain state via TMS. In particular, we 
consider the problem of driving the brain from a resting state 
to one associated with high-level cognitive operations. 

Brain resting states are associated with high activity in 
the posterior cingulate (PC) and inferior parietal (IP) cortices 
whereas active states are associated with high activity in 
the rostral middle frontal (RMF) and superior parietal (SP) 
cortices p4) , p5| . In Figure]^ we present the initial y^ and 
target yt signals, where the activity corresponding to the eight 
regions mentioned - left and right versions of each cortex - 
is highlighted with larger markers. In order to drive the brain 
from Yi to y* we consider a MN-MT seeding scheme with six 
seeding nodes. Since it is unclear how to implement a low- 
pass Alter in a human brain, we consider that each seeding 
node injects eleven values, totalizing P = 66 seeding values 
permitting the recovery of the target signal after the seeding 
phase without the need of a posterior Altering phase. Notice 
that throughout the paper we assumed the initial signal y^ to 
be zero, meaning that there is no signal present on the graph 
before the reconstruction process. However, our model can 
accommodate for y^ different from zero. To see this, if the 
seeding phase lasts t instants, then we can design our seeding 
values to recover the signal yr = y* — S^^^y^ in the original 
formulation so that the negative term cancels the effect of the 
seeding phase on yt and the target signal yt is recovered. 

We consider n oisy in jections following the constant SNR 
model in Section VTB for a = 10“^. Denoting by z the re¬ 


constructed signal, deflne the reconstruction error as e = ||z — 
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ytib/llyrib- We compute e for every possible combination 
of seeding nodes. Given that the seeding values are induced 
by TMS, we discard as possible seeding nodes the regions 
inaccessible by TMS like the ones located in the medial cortex 
and subcortical structures. After discarding inaccessible ROIs, 
the six seeding nodes can be chosen out of 38 possible ROIs, 
amounting to 2,760,681 possible configurations. In Figure [6a| 
we present a histogram of the reconstruction error for different 
seeding configurations where we only show those attaining 
errors below 0.1. The red bar in this histogram corresponds to 
the 1,611 configurations that achieve the lowest reconstruction 
errors. In Figure 6b we present the frequency of appearance 
of each ROI in these 1,611 robust seeding configurations. 
The regions with zero appearances correspond to the ROIs 
inaccessible to TMS, however, among the accessible regions 
the frequency of appearance is not uniform. For example, the 
left inferior parietal cortex in position 41 appears 583 times 
whereas the right bank of the frontal pole in position 1 is 
only used 50 times. In Figure 6c we depict the 6 regions 
more commonly used in robust seeding configurations. Notice 
that both the left and right versions of the Pars Orbitalis are 
commonly used as seeding nodes, suggesting the importance 
of this region for robust brain state induction. 


VIII. Conclusions 

A novel approach for the recovery of bandlimited graph 
signals - that admit a sparse representation in the frequency 
domain - was proposed. The focus was not on estimating 
an unknown graph signal but rather on inducing a known 
bandlimited signal through minimal actions on the graph. 
These actions referred to signal injections at different seeding 
nodes, which then percolate through the graph via local 
interactions described by a graph filter. Restrictions on the 
number of seeding nodes and the amount of injections at each 
node gave rise to three different reconstruction schemes and 
their performance in noiseless and noisy settings was analyzed. 
For the noiseless case, we showed that a AT-bandlimited signal 
can be recovered using K injections followed by a low- 
pass filter in the (graph) frequency domain. In contrast to 
classical time-varying signals, it was also shown that if the 
seeding nodes inject the values of the original signal in those 
nodes, perfect recovery is not feasible. For scenarios leading to 
imperfect reconstruction, we analyzed robust seeding strategies 
to minimize distortion. Finally, the different reconstruction 
schemes were illustrated through numerical experiments in 
both synthetic and real-world graph signals. 
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